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DTIO© 


•  The  nonlinear  thermoelastic  equations  are  examined, 
considering  large  deformations  and  temperature-dependent 
material  parameters.  Emphasis  is  on  high-speed  sliding 
contacts,  wherein  high  temperatures  localize  near  the 
contact  surface,  contributing  to  significance  there  of 
coupling  terms  neglected  m  elementary  theory. 

r 

-CMEtCLATURE  1 


Coefficient  matrices  for  governing  equations 
Combinations  of  elastic  parameters,  Eqs. (20); 
alternatively,  complex  coefficients  defined 
following  Eqs. (10),  (k •  1,  ...,  5) 

Complex  coefficients  of  material  inhomogene¬ 
ity,  Eqs.  (10),  n»  1,2 
Half  contact  width,  m 

Vector,  nonhomogeneous  part  of  difference 
equations 

Boussinesq-Papkovitch  functions,  j«  0,1,2 

Specific  heat  at  constant  deformation  (or 
at  constant  volume  for  subscript  v),  J/(kg'K) 
Value  of  cg  at  reference  temperature  tQ 

Lagrangian  strain  components 
Young's  modulus.  Pa 

Thermal  dependence  of  cE,  cE/cQi  alternative¬ 
ly,  scalar  function  of  v 
Thermal  dependence  of  k,  k/kQ 

Gradient  of  f(U),  Sf/3U 
Hessian  matrix'’of  ff  (U) ,  Eq.  (15) 

Finite  difference  grid  spacing  for 
k  -  1, 2 

Thermal  conductivity,  w/(m*K) 

Value  of  k  at  reference  temperature  Tq 

2a,  contact  width,  m 

Mechanical  equivalent  of  heat,  1.0  N'm/J 
PecletMumber,  defined  following  Eq. (7) 

Traction  coefficients,  Eqs. (20),  pa 
Contact  pressure,  Pa 

Heat 'flux,  lagrangian  description,  W 
Nonhomc jeneous  part  of  complex  traction 
boundary  condition,  Eqs. (10) ;  alternatively, 
quadratic  scalar  function  of  U,  Eq. (18) 
Temperature,  °C 

Reference  temperature,  293  K  (528°R),  unless 
indicated  otherwise 
Oimensionless  temperature,  (T-T0)/T 
Time,  s 

Solution  vector 

Oimensionless  displacement  components  u^/i 
Displacement  components,  k  •  1, 2 

Speed  of  contact  sliding,  m/s 
Generalized  wave  speeds,  Eqs. (19),  m/s 

Generalized  analytic  function,  w(z,z) 
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Xk  »  space  coordinates,  Lagrangian  description, m 
z  =■  Complex  variable,  +  i£ ^ 

y  >  Coefficient  of  linear  thermal  expansion,  1/K 
i  •  Coefficient  of  bulk  thermal  expansion  ■ 

(3\  ♦  2u)0r,  Pa/K 

•  Dilatation  rate,  conventional  (small  strain) 
notation,  1/s 

•  Transformed  space  variables,  k  •  1, 2 

7  a  Boundary  layer  (stretched)  thickness 
coordinate 

>.,u  ■  Lame  parameters,  Pa 

■  Friction  coefficient 

v  »  Poisson's  ratio 

v^  •  Third  order  elastic  parameters,  notation  of 
Toupin  and  Bernstein,  k«  1,2,3 

•  Dimensionless  space  variables,  k^/1,  k  •  1, 2 

•  Reference  density,  kg/m3 

£.  •  Components  of  Piola-Kirchhof f  Stress, 

1J  First  Kind,  Pa 
t,y  •  Analytic  functions 

i p,ip  >  Generalized  analytic  functions 

«  Generalized  Cauchy  kernels,  Eq.  (11) 

•  Exponential  factor  for 

»  Relaxation  factors,  Eq. (14) 

INTRODUCTION 

Generation  of  large  amounts  of  heat  near  the  sur¬ 
face  of  a  thermoelastic  body  subjected  to  fast-moving 
contact  loads  is  a  generally  recognized  occurrence. 
Important  in  such  a  situation  is  the  thermally-induced 
deformation,  affecting  the  contact  stress  distribution 
which,  in  turn,  affects  the  heat  input.  Stabilization 
by  limitation  of  temperatures  and  deformations  to 
steady-state  values  follows  only  if  deformations  are 
balanced  by  wear  [1,2]  leading  to  maximum  possible  con¬ 
tact  after  "wearing-in." 

Simple  contact  configurations  permit  temperature 
measurements  by  either  noncontacting  infrared  devices 
or  thermocouples  embedded  without  disturbing  temperature 
flow  (1-3);  surface  deformations  are  inferred  from 
laser/photocell  system  measurements  of  surface  clear¬ 
ances  (3) .  Such  measurements,  however,  are  insufficient 
for  complete  descriptions  of  temperatures  and  deforma¬ 
tions  and  must  be  supplemented  by  analysis.  Accurate 
analytical  predictions  aid  interpretation  of  experi¬ 
mental  results  and  provide  guidance  for  situations 
wherein  measurements  are  not  feasible. 

Conventional  thermoelastic  analyses  (linear,  un¬ 
coupled  thermal  and  elastic  effects)  often  give  tempera¬ 
ture  predictions  at  variance  with  results  of  wear  tests. 
Under  examination  herein  is  the  significance  of  thermo- 
elastic  nonlinearities,  viewed  toward  explaining  such 
discrepancies.  The  following  observations  arc  pertinent 
in  this  respect: 
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X.  Apparent  temperature  jumps  across  tha  contact 
surfaca  and  othar  a noma lias  ara  predictad  in  extrapola¬ 
tions,  by  conventional  theory,  from  near-surface  em¬ 
bedded  thermocouple  readings  [3] , 

2.  Neglect  of  theraoelastic  coupling  is  usually 
based  on  considerations  that  tacitly  assume  only  moder¬ 
ate  heat  generation  (4).  Temperatures  of  several 
hundred  degrees  generated  in  high-speed  sliding,  how¬ 
ever,  can  render  the  coupling  terms  significant. 

3.  Usual  treatments  of  thermoelastic  coupling  con¬ 
sider  only  dilatational  effects.  Significant  tempera¬ 
ture  buildups  have  been  generated  experimentally,  how 
ever,  under  repeated  cycling  of  elastic  deviatoric 
strains  [S] . 

4.  Variation  of  thermal  parameters  with  temperature 
can  significantly  affect  heat  conduction  in  sliding  con¬ 
tacts,  as  analysis  demonstrates  (6]. 

The  reality  of  irreversible,  often  large,  deforma¬ 
tions  and/or  strains  accompanied  by  heat  generation 
under  inelastic  conditions  lad  long  ago  [7]  to  recogni¬ 
tion  of  the  importance  therein  of  thermomechanical 
coupling.  In  thermoelasticity,  however,  a  typical  argu¬ 
ment  [4]  considers  only  the  simplest  coupling,  involving 
a  dilatational  term  in  the  heat  conduction  equation  with 
magnitude  governed  by  the  factor  3  (3o7g/fjc, )  («„,  /3aT) . 
Using  isothermal  values  of  3,  a  c,  for  a  typical  steel 
and  for  aluminum  and  366  K  (660°R)  for  Tq,  the  factor 
3  <5oT0/pc, )  is  evaluated  as  0.01793  and  0.03714, 
respectively;  it  is  concluded  that  coupling  is  negli¬ 
gible  if  (etk/3oT)«20  and  this  is  assumed  to  be  the 
case.  On  the  other  hand,  high-temperature  values  [8-10] 
of  3,  or,  c,  for  a  typical  steel  at  565°C  (1230°  P)  -  an 
attainable  dry  sliding  surface  temperature  [6]  -  give 
0.1576  for  3  (SorT^/pc, ) .  Thus,  even  if  (t|,,/3a*)  is  only 
of  order  1.0  the  coupling  terra  in  this  case  has  an  order 
-  16%  of  the  basic  term  gc.T.  The  apparent  importance 
of  this  simplest  coupling  terra  suggests  examining  the 
significance  of  additional  relevant  coupling  terms,  as 
well. 

THEORY 


Figure  1  Schematic  Diagram  of  the  Half-Space 
Subjected  to  Mechanical  and  Thermal 
Loading 

frictional  force  generates  dissipative  energy  that  pro¬ 
duces  a  heat  source  with  flux  Q(Xl)  »  Vp,  p(Xx  )/M;  the 
remainder  of  the  surface  is  assumed  insulated.  Pres¬ 
sure  p(Xx )  is  considered  to  be  the  force  on  the  deformed 
surface  per  unit  area  of  undeformed  surface,  since  the 
deformed  shape  is  not  known  a  priori.  Thus,  the  applied 
loads  are  surface  tractions  related  to  piola-Kirchhoff 
stresses  of  the  First  Kind  [13-15],  Similarly,  heat 
flux  Q (Xx  )  is  that  referred  to  the  undeformed  solid 
[13-15] . 


Analytical  Model 

Evidence  that  high-speed  sliding  contact  situations 
produce  temperatures  and  material  degradation  localized 
mainly  near  the  contact  surface  permits  consideration  of 
a  semi-infinite  solid  with  contact  loads  moving  along  a 
plana  boundary.  Applying  the  plane  boundary  idealiza¬ 
tion  to  problems  involving  curved  surfaces  gives  an 
error  in  the  governing  equation*  of  order  PguaR“l  ,  with 
PQ  the  Pec  let  Number  and  R  the  ratio  of  tha  surface 
radius  of  curvature  to  a  typical  length  (e.g.,  the  sur¬ 
face  contact  width) .  In  addition,  only  a  two- 
dimensional  model  is  needed  because  high-speed  sliding 
leads,  essentially,  to  suppression  of  heat  flow  in  the 
surface  direction  perpendicular  to  sliding  [11] .  The 
model  utilized  is  a  plane  strain  one,  applicable  to  a 
large  contact  lingth  in  tha  latter  direction. 

Figure  1  shows  the  undeformed  solid,  subjected  to 
mechanical  loads  p(Xi)  and  u«  P  )  over  the  surface 
region  0sx1  <  2a,  X,  ■  0,  |x,  |  <•,  moving  with  constant 
speed  V  relative  to  the  bulk  solid;  the  latter  load  is 
a  frictional  force  per  unit  of  material  thickness  in  the 
X, -direction.  Coordinate  system  (k  •  1,2,3)  is  fixed 
with  respect  to  the  material,  whareas  system  x*  is  fixed 
with  respect  to  the  leading.  The  saterial  is  considered 
to  be  isotropic,  free  of  body  forces,  and  to  have 
temperature-dependant  material  parameters.  In  addition, 
the  friction  coefficient  is  generally  a  nonlinear  func¬ 
tion  of  velocity,  linear  asymptotically  at  high  sliding 
speeds.  Considering  steady-state,  high-speed  sliding 
fixes  the  friction  coefficient  and  rules  out  self- 
sustained  ('chattering*)  oscillations  [12],  The  surface 


Governing  Equations 

With  respect  to  reference  frame  X*,  convected  with 
the  surface  loading  according  to  X,  ■  ^  +5klvt,  dis¬ 
placements  and  temperatures  are  quasi-stationery  [6,11], 
Thus,  e.g.,  u*  (Xj,  Xj',  t)«uk(Xl-Vt,  Xa),  so  that 


etc.  Hence,  velocities  and  accelerations  referred  to 
material  coordinates  can  be  expressed  entirely  in 
terms  of  convected  coordinate  Xx,  without  explicit 
dependence  on  t.  Application  to  tha  nonlinear  thermo¬ 
elastic  equations  referred  to  material  coordinates  xj 
[15]  permits  their  expression  in  terms  of  the  convected 
coordinates.  The  two  momentum  equations  and  the  heat 
equation  (after  rearranging  and  nondiswnsionalising) 
are  then,  respectively. 


Similarly,  the  boundary  conditions  take  the  form: 


C22(^l’0)  “  [PNllUl,  1  +  PN22U2,2  * PN12U1, 2  +  PN21U2, 1 
-ST']  =.-f*P<5l\  °S5lSl\ 

@  s2*o  1  0,5^  to,  u  ; ' 
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E12(*l'0)  *  [PS12U1,2  +  PS21U2,  l!52-0 
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The  coefficient  terms  v^, ,  p,,,  appearing  in  Eqs. (2)  - 
(6)  are  defined  in  Appendix  A.  In  addition,  dimension¬ 
less  variables  are  U,  »  u ,/jt,  5,  «X,/Z,  T'  »  (T  -  Tj  )/%  , 
f  (T'l  «c,  (T')/q>,  g(T' )  «  k(T')/K] ,  and  abbreviations 
(  ),,«&<  )/3£, ,  (  ) '  -  a (  )/3T',  i«2a  apply.  P0  is  the 
Peclet  Number  VfpoCg/%.  The  v?j,  all  have  dimensions 
of  (velocity)3  and  vlu,  viaa  generalize  the  classical 
isothermal  speeds  of  dilatational  and  shear  waves, 
respectively. 

Equations  (2)  -  (4)  differ  from  versions  given  pre¬ 
viously  [IS]  by  allowing  nonlinear  (rather  than  simply 
linear)  thermal  variations  of  material  parameters  and  by 
not  neglecting  T'  compared  to  1  in  [v^n,  ...,  ^).  How¬ 
ever,  dependence  of  thermal  conductivity  on  strain  (IS), 
apparently  never  quantified  experimentally,  is  neglected 
herein;  further,  the  formulation  allows  for  thermal  de¬ 
pendence  of  third  order  elastic  parametsrs  vj,  v,,  \jj , 
but  lack  of  data  forced  use  of  only  isothermal  values 
in  calculations. 


Boundary  layer 

Typical  sliding  contact  speeds  are  much  mealier 
than  elastic  wave  speeds  (e.g.,  by  orders  *  10*  vs. 

10s  in. /see  for  a  typical  steel),  although  they  can 
yield  large  values  of  P0.  Thus,  terms  tl  -  (V/vn  j  )a  Jv?;l , 
etc.,  in  momentum  Iqs. (2),  (3)  can  be  considered  to  be 
of  order  0(1).  Second  derivatives  of  T'  in  the  heat 
Eq. (4),  on  the  other  hand,  are  multiplied  by  Fgl  which 
becomes  mull  for  large  >b .  Hence,  the  term  containing 
these  derivatives  becomes  vanishingly  small  for  large  P0, 


except  in  regions  where  (T.'n  ♦T'aa)  is  O(P0).  A  first 
approximation  to  the  solution  calculated  from  this  re¬ 
duced  differential  order  problem  would  be  incapable  of 
satisfying  all  specified  thermal  boundary  conditions  - 
in  particular,  Eq.  (7)  for  £j  •  0  vs.  T'-O  as  £,-•. 
Analysis  for  the  relevant  characteristic  curves  clari¬ 
fies  this  difficulty.  The  direction  cosines  n,,  n,  of 
such  a  curve,  related  by  rg  »  n;>  -  1,  are  found  [16]  from 
the  relation  det  (A) *  0,  where  A  is  the  coefficient 
matrix  for  the  vector  [UjUbT')  (transposed)  in  Eqs.  (2)  - 
(4),  with  d/d£,  replaced  by  n,  ( i  -  1,2)  [17).  Th*  re¬ 
sulting  algebraic  equation  is  of  fifth  order  (rather 
than  sixth  for  the  full  problem)  and  has  a  real  root 
•  0,  giving  curves  £3  •  constant  generated  by  reduction 
in  order  of  the  heat  equation.  Vanishing  temperatures 
prescribed  (as  Cauchy  data)  as  t  -  *  should  then  define 
temperatures  throughout  the  entire  semi-infinite  domain 
considered.  The  discrepancy  so  generated  as  £3  -0  indi¬ 
cates  the  existence  there  of  a  boundary  layer  with 
5T'/S£3  large. 

The  remaining  quartic  polynomial  generated  by 
det (A)  » 0  is  presumed  to  not  yield  real  roots  (real 
characteristics)  associated  with  elastic  waves  because 
of  the  smallness  of  V/v, ,,  noted  above.  Another  possi¬ 
bility,  however,  is  that  the  momentum  equations  may  fail 
to  be  elliptic  for  displacement  gradient  magnitudes 
exceeding  certain  bounds,  as  discussed  recently  [18, 19) 
for  isothermal  plane  finite  elastostatics;  treatment  of 
thermal  effects  herein  gives  revised  criteria.  Consid¬ 
eration  in  what  follows  is  restricted  to  the  assumption, 
considered  plausible  for  typical  sliding  situations, 
that  displacement  gradients  and  temperatures  are  only 
moderately  large  away  from  the  contact;  thus,  the  lattar 
possibility  fails  to  materialize.  Then  it  is  inferred 
(17)  that  the  quartic  has  only  imaginary  roots,  so  only 
the  above-mentioned  boundary  layer  near  £3-0  occurs. 

It  can  be  shown  that  the  boundary  layer  identified 
above  has  a  thickness  of  order  0(Pjl/i).  correspond— ;g 
to  stretched  thickness  coordinate  I]' Eg  £3.  Thermo- 
elastic  behavior  therein  is  controlled  mainly  by  reduc¬ 
tion  of  the  heat  equation  from  elliptic  to  parabolic 
type  (having  characteristics  perpendicular  to  the  con¬ 
tact  surface  rather  than  parallel  to  it  as  in  the 
reduced  case  indicated  above)  and  of  the  momentum  aqua¬ 
tions  to  ordinary  differential  variety.  Outside  the 
boundary  layer,  temperatures  are  generally  suppressed. 

In  contrast,  deformations  are  not  so  restricted,  but 
even  to  first  approximation  are  associated  with  a  po¬ 
tential  theory  problem  involving  elliptic  equations  (see 
below).  Significant  boundary  layer  contributions  to  the 
deformations  are  exhibited  only  in  higher  approximations. 

Detailed  analysis  [17]  shows  that  the  boundary 
layer  is  multi-structured  and  contains  subregions  smaller 
than  0(Pgl/3)  in  which  various  approximations  apply. 

For  example,  in  small  regions  of  size  0(Pgl)  about 
£t  •  0, 1  no  approximation  is  possible  and  the  full  ellip¬ 
tic  equations  (and  potential  solutions)  govern.  The 
boundary  layer  structure  is  shown  schematically  in 
Figure  2;  tome  features  are  analogous  to  those  for 
boundary  layer  flow  past  a  flat  plate  [20,21]. 

Potential  Theory  Solutions 

In  theory,  it  it  possible  to  obtain  a  solution  by 
successive  iteration  between  the  heat  and  momentum 
equations,  considering  each  in  turn  separately.  This 
concept  gives  at  least  qualitative  insight  into  the 
solution  character.  Treating  the  nonlinearities  in  the 
heat  equation  as  a  nonhomogeneous  term  P* (£i,£a>  and 
removing  thermal  dependency  of  k  (i.e. ,  of  g)  by  a 

transformation  [6]  sp  »  J  g(T")dT",  due  originally  to 
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Figure  2  Schematic  Representation  of  the  Thermal 
Boundary  Layer  Structure 


method  [25] ,  the  difference  being  that  the  present  com¬ 
plex  functions  are  displacement,  not  stress,  functions. 
Thus  the  (usual  and  generalized)  analyticity  properties 
utilized  herein  follow  from  the  (nonhomogeneous)  bihar¬ 
monic  character  of  the  Galerkin  vector,  rather  than  that 
of  an  Airy  stress  function;  complicated  considerations 
of  displacement  compatibility  are  thus  avoided  by  this 
displacement  formulation. 

The  complex  form  of  boundary  conditions  (5),  (6)  ie 

•l(Ci>«/<C1>*»2<Ci>«T<s£7  *  C1*’’Tc]T 


♦  7rT^T-a3<;1)s<5[) , 
a2<Ci)4'(Ci)*a1(C1)?T(CjT  * 


Kirchhoff  [22],  permits  expression  of  the  solution  for 
cp  as  (j  ■  <p,  +  %  ,  where 

i  }  y  <5i-5i> 


•  Ko(f  (si 

*  2n  I  6  d5i  J 

-o»  o 

P 

•  Wf  K-en>2+<V^>2)  <9> 

+  Ko(t  ^5X-5:)2>(C2-S2)2)>*(!{,5i.a?' 


and  l^j(fl)  is  the  modified  Bessel  function  of  the  second 
kind,  order  zero.  Equation  (8)  is  equivalent  to  a  pre¬ 
vious  result  [22,  p.269],  obtained  by  superposition  of 
heat  sources,  for  only  surface  heat  application  by  uni¬ 
form  heat  flux  over  a  strip  area.  Asymptotic  considera¬ 
tions  [17]  show  reduction  as  P0  -■»  of  Eqs.  (8),  (9)  to 
forms  [6]  having  kernels  (Si  -  5a  )"1  /a exp[-0] +  T|*)3/ 

4(ft  -?{)]  with  T|*  •  0  for  g>,  and  T]*-±T|'for  tp,  ;  inte¬ 
gration  for  Si  in  then  extends  to  0,  5, ,  1  for  5t  <  0, 
0  s  Si  s  1,  Si  >  lt  respectively, and  in  always  to  Sj 
(other  ranges  of  Si  yield  exponentially  small 
contributions) . 

Thermal  dependence  of  the  elastic  materiel  param¬ 
eters  can  be  mollified  by  transforming  space  variables 
in  the  momentum  equations  in  analogy  to  that  Indicated 
above.  Treating  remaining  nonlinearitias  end  the  accel¬ 
eration  terms  in  these  equations  as  equivalent  body 
forces  yields  equations  amenable  to  elastic  potential 
solution.  Thus,  a  formulation  similar  to  Papkovitch  [23] 
shows  that  the  Calerkin  vector  satisfies  a  nonhomogene¬ 
ous  biharmonic  equation  and  identifies  Bousslneaq- 
Papkovitch  functions  £(, ,  Mf  separable  into  homogen¬ 
eous  (harmonic). end  particular  components  »iH ,  B)„ 
end  Bi,,  ■», ,  respectively.  The  former  are  related 
to  analytic  functions  l(z),  T(s)  of  a  eoaplax  variable 
*  *  Ci  ♦  iCe  <Ci » it  transformed  space  variables)  by 
Bn  ♦  i^,  -2(1  +  ^)  t  (a)  and  M^/Ba-  (1  ♦  A,  )T(1),  with 
Ag  a  constant.  Relations  of  the  same  type  applied  to 
B,, ,  Bir ,  define  generalised  analytic  functions 
t,  (s,2>,  1,  (s,i)  in  the  sense  treated  by  Vekua  [24] . 

As  a  result,  the  complex  displacement  ut  ♦  iu,  can  be 
expressed  as  a  generalisation  of  the  usual  form  [25, 
p.1121,  the  additional  terms  being  due  to  t, (s,l), 
f»(s,S).  Traction  boundary  conditions  (5),  (6)  are 
treated  effectively  by  a  variant  of  Muskhelishvlli’ s 


+  Y' <cx>  •a3n1)s<c1)  ,  (10) 


with  S  (Ci )  and  its  complex  conjugate  S(J, ),  involving 
complicated  combinations  of  surface  Piola-Kirchhoff 
stresses  and  linear  and  nonlinear  displacement  gradient 
and  temperature  terms,  assumed  known  from  previous  iter¬ 
ations.  Equations  (10)  are  analogous  to  Muskhelishvili 
[25],  Eqs.  (93.3),  (93.4),  with  the  added  complication 
that  Si  (Ci),  as(Cs)  are  surface  values  of  functions  that 
can  vary  spatially  due  to  thermally-induced  material 
inhomogeneity.  However,  the  functions  wt (z,i)  - 
Si(z,z)4'(z)  and  wa  (z,2)  *  aj  (z,  z)t '  (z)  satisfy  nonhomo- 
geneous  Cauchy-Riemann  equations  BW* /Bz » \ n , 

\  -  Stlog  aj  (z,z)]/dz,  (k-1,2),  etc.  and,  hence,  are 
specializations  of  generalized  analytic  functions  [24] . 
Consequently,  the  w^  (z,z)  possess  generalized  Cauchy 
kernels  (z,Zt  ),  viz.  : 


nik(*»*B* 


•"ik^'V 


*B'* 


(k-1,2) 


lk 


(z'  -  z)  (zB-  z') 


atidCi 


(11) 


where  Zt  denotes  a  point  on  the  boundary  and  domain  D 
is  the  semi-infinite  region  under  consideration. 

Terms  Cx  S  (Ci  ) .  t '  (Ci  )  in  Eqs.  (10)  are  boundary  values 
of  functions  z4  (z), _T  (z)  holomorphlc  in  the  lower 
half-plane  and  w,  (C.  C,  )  ■  w .  C.  )  are  boundary  values 
of  functions  wl(z/z),  w3(z,z)  analytic  in  a  generalized 
sense  in  the  upper  and  lower  half-planes,  respectively. 
Application  of  the  usual  Cauchy  integral  theorem  [25] 
and  a  generalized  Cauchy  formula  [24,  p. 175]  thus  yields 
relations  for  4'(z),  f'(z); 


a^l  z,l)4'(z) 
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"  2m  J  c,  -  s  V^i1*  l^i,d^l 
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Q*({i)  is  che  dimensionless  surface  heat  flux,  given 
by  the  left-hand  side  of  iq. (7). 
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Equations  (12)  generalize  results  of  Muskhelishvili  125], 
Eqs. (93.6),  (93.7),  and  demonstrate  relevance  of  an 
elastic  potential  solution  for  nonhomogeneous  materials; 
the  first  is  a  linear  integral  equation,  theoretically 
solvable  [17] . 

NUMERICAL  ANALYSIS 
General 

The  complexity  of  Eqs. (2)  -  (7)  rules  out  their 
complete  analytical  treatment.  Thus,  their  solution 
must  be  accomplished  numerically,  using  iteration.  To 
enhance  convergence,  simultaneous  operations  with 
Eqs. (2)  -  (7,  seemed  advisable.  An  iterative  finite 
difference  scheme  for  this  purpose  idealized  the  semi- 
infinite  region  by  a  finite  rectangular  grid,  Figure  3, 
guided  by  analytical  results  mentioned  above.  Detailed 
application  of  these  results  to  determine  grid  size 
and  spacing  is  described  elsewhere  [17],  Briefly, 
however,  known  exponential  decay  of  temperature  in 
front  of  the  contact  (£.  <  0)  and  through  the  boundary 
layer  thickness  facilitates  specification  of  vanishing 
temperatures  and  of  displacements  obtained  by  linear 
analysis  at  the  boundaries  £i  »-  0.5  and  §,  «  10  PS1'2 
(T)  -  10 ) .  Behind  the  contact  (5l>0),  temperature  decay 
is  merely  algebraic,  so  choosing  ^  large  enough  for 
specification  of  (approximately)  vanishing  temperatures 
there  is  not  feasible. 

If  Eqs. (2)  -  (4)  were  of  either  parabolic  or 
parabolic-hyperbolic  character,  formulation  as  an 


initial  value  problem  in  f ,  utilizing  forward  integra¬ 
tion  through  explicit  finite  differences,  would  elimi¬ 
nate  the  need  for  boundary  values  at  ,,,  .  The  basic¬ 
ally  elliptic  character  of  Eqs. (2),  (3)  and  the  ellip- 
ticity  of  Eq.  (4)  about  “0,1,  Cs  »  0,  together  with 
stability  questions  inherent  in  explicit  schemes,  how¬ 
ever,  indicated  treatment  as  purely  a  boundary  value 
problem,  using  an  implicit  difference  scheme.  An  effort 
to  limit  the  number  of  grid  points,  dictating  a  compro¬ 
mise  between  high  grid  detail  and  the  magnitude  of 
Si,.,  used,  resulted  in  £la>,  «  2. 5  and  specification 
there  of  temperatures  and  (approximately)  vanishing 
displacement  gradients  Uy  f,  u2  ,  obtained  by  linear 
analysis.  Boundary  conditions  '(5)  -  (7)  utilized  a 
Hertzian  distribution  of  contact  pressure  p(£, ) .  The 
elementary  solution  for  temperatures  then  involves  con¬ 
fluent  forms  of  a  hypergeometric  function  of  two  variables 
[17],  but  the  usual  result  [22]  based  on  a  uniform 
average  p(^)  over  the  contact  gives  an  asymptotically 
accurate  approximation  at  >2.5. 

Finite  Difference  Form  of  the  Governing  Equations 

Development  of  the  finite  difference  approximation 
of  the  governing  equations  followed  the  form  of  these 
equations  exhibiting  isolated  highest  derivatives  of 
the  dependent  variables.  Thus,  difference  operators 
were  applied  directly  to  Eqs. (2)  -  (7).  Use  of  nonuni¬ 
form  spacings,  Figure  3,  however,  necessitated  more 
general  forms  of  difference  operators  than  are  commonly 
available.  Thus,  approximation  of  derivatives  with 
error,  order  of  the  grid  spacing  squared,  generally 
involved  12  grid  points,  whereas  uniform  spacing 
requires  nine. 

The  3x  34  X  10“  1020  difference  equations  have  the 
matrix  form 

A (U)U  •  B  ;  (13, 

complicated  construction  of  Eq. (13)  is  detailed  else¬ 
where  [17] .  Equation  (13)  is  nonlinear  because 
matrix  A  is  a  function  of  solution  vector  U,  consisting 
of  Uj,  Hj,  and  T'  taken  successively  at  each  point  (J,J> 
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Figure  3  Sample  Finite  Difference  Grid  for  P„  ■  1000,  Illustrating  variable  Grid  Spacing  Needed 
Through  the  Boundary  Layer  Thickness  and  at  the  Contact  Region  Leading  and  Trailing  Edj 
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left  Co  eight  in  row*  1  through  10.  Matrix  A  is  unsynr- 
metric  but  has  banded  character,  with  41  diagonals  on 
which  nonzero  elements  can  occur.  Thus,  A  has  (1020)3  • 
1,040,400  elements,  but  only  27,216  are  nonzero.  Gen* 
eration  of  the  difference  grid,  calculation  of  coeffi¬ 
cients  for  the  difference  equations,  construction  of 
matrix  A  and  vector  B  and  solution  of  Eq.  (13)  were 
programmed  for  IBM  3033  computation.  NAG  (26]  sub¬ 
routines  F03AJF  and  F04APF,  designed  for  large,  sparse 
matrices,  were  used  to  triangularize  matrix  A  and  solve 
Eq.  (13)  by  Gaussian  elimination.  Economy  of  storage 
space  required  identification  of  only  nonzero  elements 
of  A  and  their  locations  in  A.  However,  nonzero  ele¬ 
ments  occur  in  four  diagonal  bands  11,  14,  11,  and  5 
elements  in  width,  separated  by  diagonal  null  bands  each 
102  zeros  in  width.  Hence,  storage  requirements  for 
the  triangular  factors  considerably  exceed  that  for  the 
original  27,216  nonzero  elements;  the  actual  number  of 
elements  generated  during  triangularization  is  about 
100,000. 

Numerical  Calculations 

Numerical  calculations  treated  a  typical  contact 
situation,  based  on  Figures  1  and  3,  wherein  values  of 
sliding  speed  (0.8611  m/s),  Peclet  Number  <P0«  1000), 
average  heat  flux  (4.43X1Q7  N/(m*s))  and  thermal 
parameters  were  chosen  to  correspond  to  those  of  an 
earlier  study.  That  study  [6]  considered  temperature- 
dependent  thermal  parameters;  thus,  an  envisioned  com¬ 
parison  of  results  obtained  therein  for  SAE  1010  steel 
with  those  sought  in  the  present  study  for  the  same 
material  offered  to  clarify  the  effect  of  the  added 
nonlinaarities.  Quadratic  variations  of  1,  u  with  T' 
were  used,  based  on  tabulated  information  (101  for  E, 
v ;  similar  variations  were  employed  for  vl(  va,  u, , 
but  only  isothermal  values  for  iron  (27]  could  be 
located  as  the  closest  applicable.  Although  ur  is  a 
complicated  function  of  V,  T'  and  other  factors,  in¬ 
formation  from  several  sources  (17)  suggested  a  tem¬ 
perature  dependency  of  error  function  type  (S-jhaped) 
with  y,  ranging  from  about  0.115  for  20  -—  205  C  to 
about  0.400  above  -  595  C.  Hertzian  pressures  used 
corresponded  to  a  contact  load  of  5.838xl(f  N/(ra-s). 

Direct  searches  were  made  for  values  of  U, ,  U, ,  t' 
at  the  340  grid  points.  Ideally,  these  would  lead  to 
solution  of  nonlinear  Eq.  (13)  by  a  sequence  of  Newton- 
Raphson  iterations  (Newton's  Method  in  n  dimensions, 

•  1020  for  3x  340  unknowns),  corresponding  to  succes¬ 
sive  linear  solutions  with  updating  of  matrix  coeffi¬ 
cients.  In  theory,  this  process  is  repeated  until 
convergence  is  attained,  but  modification  utilizing 
underrelaxaeion  [28]  was  found  necessary  to  enhance, 
or  possibly  guarantee,  convergence.  Thus,  for  itera¬ 
tion  k  an  accepted  value  of  each  component  U,  of  U  isi 

Uilt>  “*i®i°  ♦  U-«t)o‘k'l>  ,  i-  1,  ...  ,  1020  .  (14) 

In  Eq.  (14),  the  (?,*’  are  components  of  the  current 
solution  vector  and  the  (/,  are  corresponding  com¬ 
ponents  from  the  previous  iteration.  The  u>,  are  relax¬ 
ation  factors,  assigned  uniform,  distinct  values  iu, , 
is,,  tie  for  0ll  Uj,  T',  respectively,  variable  with  k, 
although  variation  of  the  s,  with  location  could  have 
been  allowed  also.  Underrelaxation  rather  than  over- 
relaxation,  raters  to  choices  w,  s  1.  dictated  by  solu¬ 
tion  character;  thus,  Eq. (14)  indicates  starting  values 
obtained  by  interpolation  between  the  results  of  the 
previous  two  iterations. 

The  above-mentioned  searches  aimed  at  establishing 
the  general  range  of  U,,  li,,  T'  wherein  the  true  solu¬ 
tion  U  might  lie,  subsequent  iteration  being  deferred 
until*‘convergencs  in  wall  steps  seamed  reasonably 


certain.  Progress  was  limited,  however,  because  large 
storage  requirements  for  triangularization  of  A  made 
each  linear  solution  expensive;  moreover,  setting  relax¬ 
ation  factors  is,  by  trial  and  error  proved  inefficient. 
Thus,  trial  starting  values  of  the  U,  generated  by 
Eq.  (14),  if  not  close  to  those  of  the  true  solution, 
defined  "long"  Newton-Raphson  tangents  (hyperplanes  in 
1020-degroe  Euclidian  space),  giving  new  values  of  the 
U.  from  Eq. (13)  often  far  from  the  starting  values. 
Apparently,  convergence  speed  was  dependent  on  grid 
point  location,  seemingly  analogous  to  locally  varying 
degree  of  stability  in  forward  integration  of  nonlinear 
parabolic  equations  [28,  p.325].  Convergence  speed  was 
undoubtedly  dependent  on  local  grid  spacing  ratio  l\/hj, 
also;  although  this  ratio  is,  by  far,  not  as  critical 
in  boundary  value  problems  as  in  initial  value  formula¬ 
tions,  large  values  of  h,/h,  gave  locally  anomalous 
results  (mainly  for  T*  central  to  the  contact  surface). 
The  searches  made  indicated  that  a  particular  choice  of 
the  is,  offered  possible  convergence  at  soma  locations 
but  divergence  at  others.  Truly  effective  relaxation 
must,  then,  allow  for  distinct  s,,  m,,  iu,  at  each  of 
340  points,  but  setting  all  such  values  by  trial  and 
error  would  be  impossible. 

A  more  efficient  search  technique  considers  the 
quantity  A(U)U-B  of  Eq.  (13)  to  represent  the  gradient 
g  (U>  of  some  scalar  function  f(U).  Thus,  satisfaction 
of  Eq. (13)  is  equivalent  to  finding  the  global  minimum 
of  hypersurface  f(U)  *0.  A  quadratic  expression  for 
ftU)  applies  when  A  is  a  symmetric  matrix  of  constants; 
a  more  general  expression  [171  accounts  for  unsymmetry 
of  A  and  dependence  of  A  on  U.  Minimization  of  f(U) 
would  be  difficult.  Indeed,  starting  values  for  some 
searches  led  to  det  [A(U>]  <  0,  indicating  A(U)  not 
everywhere  positive  definite  in  the  space  of  U  and  f(U) 
not  everywhere  convex.  Physical  realism,  however, 
requires  restriction  co  U  giving  positive  definite  A[U). 
Synsnetry  of  A(U)  would  guarantee  n  «-1020  mutually 
orthogonal  eigenvectors,  "conjugate"  with  respect  to 
A (U) ,  spanning  the  space  of  U  and  providing  n  ■  1020 
mutually  orthogonal  search  directions  for  conjugate 
gradient  minimization'  .  Definiteness  of  A(U)  ensures 
that  none  of  the  search  directions  is  in  che  null  space 
of  A (U)  and  positive  definiteness  permits  quadratic 
termination  (theoretically,  in  n  steps)  of  the  minimiza¬ 
tion  algorithm  [29] . 

Difficulties  related  to  unsymmetry  of  A  can  be 
avoided  by  utilizing  the  Hessian  matrix 
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Practical  calculation  of  oA/au  for  Eq. (15)  utilizes 
linear  approximation  for  a  small  variation  AUi 
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For  independent  (scalar)  variations  AU,  of  individual 
U,,  iaA/au)AU  -  OA/au, ) AU,  (no  sum);  thus  Eq.  (16)  gives 
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Conjugacy  in  a  wider  sense  for  a  square,  positive 
definite  matrix  requires  only  that  the  search  directions 
be  linearly  independent  [29,  p. 153] ;  these  directions 
may  then  be  orthogonalized  by  the  Gram-Schmidt  process 
(29,  p. 162] 


Th«  cam  H(U)  is  syraietric,  indicating  possible  effect* 
iveness  of  minimising  the  quadratic  form  (with  C  a 
scalar  constant) : 

S(U)  -  t  U  •  HCUUJ- B  •  U  +  C  (18) 

where  S(U)  gives  a  local  quadratic,  convex  fit  to  f(U>. 
Hence,  if  a  search  has  located  starting  values  reason¬ 
ably  close  to  the  true  solution,  locating  the  global 
minimum1  of  S(U)  approximates  that  of  f(U).  An 
available  algorithm  [29,  pp,16S-167]  for  minimizing 
SlU)  needs  further  development  for  a  sparse  matrix 
capability  and  for  programming  in  FORTRAN.  Algorithms 
for  storing,  transposing  and  multiplying  sparse  matrices 
[30,  Chap. 3)  may  be  useful  also.  Fully  operational 
conjugate  gradient  subroutines  available  to  the  author 
[26,31-33]  either  lacked  a  large,  sparse  matrix  capabil¬ 
ity  or  required  matrix  symmetry  and  positive 
definiteness. 

Although  full  solution  by  direct  searches  was  not 
feasible,  the  results  obtained  appear  to  support  a  sig¬ 
nificant  influence  of  nonlinear  thermoelastic  effects 
in  high  temperature  situations  such  as  that  considered. 
Thus,  the  results,  although  not  fully  conclusive,  im¬ 
plied  a  maximum  surface  temperature  (at  a  0.8)  of 
roughly  780SC  (T'-2.6),  contrasting  with  about  535SC 
considering  nonlinear  thermal  effects  along  [16]  and 
about  615°  C  from  elementary  theory  [6],  Similarly,  a 
maximum  surface  outward  normal  displacement  of  roughly 
0.003  in.  was  implied.  (Corresponding  strains  as  large 
as  roughly  0.01  in. /in.  indicated  some  possible  yield¬ 
ing  for  the  steel  considered.)  Equations  (18),  (19) 
show  that  displacement  gradient  effects  begin  to  domi¬ 
nate  contributions  from  Xlti  when  |a4U.  ll,  Ia^U*  a|, 

lAiua,  ali,  l*sui.  I  *0(X*  Jii),  |a !  a,  ua ,  a  ]  - 
ott*ii),  IAiUi.J,  IajUj^I  «0(u),  etc.  For  the  steel 
considered  these  criteria  reduce  to  |ol(  ! |  -  0(0.006), 
|ua,a|  -0(0.02),  although  lack  of  high-temperature  data 
for  Vj,  va,  u,  permits  comparisons  for  only  isothermal 
parameters.  Occurrence  of  locally  negative  squares  of 
wave  speeds,  Eqs.  (8),  may  lead  to  failure  of  matrix  A 
to  be  positive  definite,  a  condition  presumably  related 
to  det(A)  <  0  noted  for  some  searches. 
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APPENDIX  -  COEFFICIENT  TERMS  FOR  EQS.  (2)  -  (6) 


PNll“ 
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